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The energy levels of hydrogen and helium atoms in strong magnetic fields are calculated in this 
study. The current work contains estimates of the ground and first few excited states of these 
systems that are improvements upon previous estimates. The methodology involves computing 
the eigenvalues and eigenvectors of the generalized two-dimensional Hartree-Fock partial differential 
equations for these one- and two-electron systems in a self-consistent manner. The method described 
herein is applicable to calculations of atomic structure in magnetic fields of arbitrary strength as 
it exploits the natural symmetries of the problem without assumptions of any basis functions for 
expressing the wave functions of the electrons or the commonly employed adiabatic approximation. 
The method is found to be readily extendable to systems with more than two electrons. 



I. INTRODUCTION 

The motivation to study atoms in magnetic fields of 
strength beyond the perturbative regime was in a large 
part due to the discovery of such fields being present in 
white dwarf stars [TJ [2J [3] and neutron stars |H \5 . The 
most commonly observed neutron stars - pulsars, have 
been observed to have magnetic fields on the order of 10 11 

- 10 13 G [BJ. Magnetars [7], which are strongly magnetized 
neutron stars, can have magnetic field strengths well in 
excess of 10 13 G. White dwarf stars on the other hand 
have somewhat less extreme fields, albeit still high, ~ 10 6 

- 10 8 G [BJ. At the high field strengths observed in these 
compact objects, the electron cyclotron energy of an 
atom becomes greater than the corresponding Coulomb 
potential energy, i.e. Hlub > Ze 2 /r [BJ. Here, lob is the 
cyclotron frequency. Thus, a Zeeman-type perturbative 
treatment of the field [8] is not possible. The structure 
of atoms in both instances however, is considerably al- 
tered from the low field case. At high field strengths, 
spherical symmetry of the atom is broken and the atom 
is stretched along the field, however, azimuthal symme- 
try remains intact and thus it has been observed that 
it is more convenient to model the atom in cylindrical 
coordinates [BJ. 

Since the 1970's this problem has been tackled by 
various researchers using different methods. Hitherto, 
the most tractable and accurate approaches have re- 
lied upon an assumed basis of functions for express- 
ing the electron wave functions of an atom in a strong 
magnetic field. Using such an assumption the prob- 
lem of the hydrogen atom in a strong magnetic field 
was tackled using either a variational approach [5] or 
by attempting to solve the Schrodinger equation directly 
[TU1 EQ ED ED HI ESI E3- Initial attempts for estimat- 
ing the energies and wave functions of different electronic 
states of the helium atom were based upon a purely vari- 
ational approach [HlllHlEaiMlEllEaESllSlES] or 
Z— dependent perturbation theory [3BJ. The majority of 
these studies however were analytical and were therefore 
limited in their applicability to a problem that was inher- 
ently more tractable numerically. The most accurate and 
reliable solutions to the helium atom in a strong magnetic 



field, thus far, involved using the Hartree-Fock (HF) tech- 
nique |27j . Employing the HF technique researchers in 
the past have calculated energy levels and wave functions 
of the ground and first few excited states of the helium 
atom and helium- like species [28j ESJ [30] . These trea- 
tises employed Landau orbitals [8] to describe the mo- 
tion of the electron perpendicular to the field. The elec- 
trons were required to reside in the ground Landau state, 
thereby simplifying calculations somewhat by restricting 
the wave functions to the so called adiabatic approxima- 
tion [25] . Ivanov [3T] in 1994 obtained similar results 
for the helium atom in strong magnetic fields using an 
unrestricted HF technique. Elsewhere, Quantum-Monte- 
Carlo methods were employed successfully for determin- 
ing the ground and first few excited states of the helium 
atom in low to strong magnetic fields [32 [331 EH]- In 
the treatises described above [28 | f29 | f30 | I3T1 . [32| [33l [34] . 
usually an approximation was employed for calculating 
both the direct and the exchange interactions between 
the electrons. Such approximations generally involved 
finding appropriate expansions that mimicked the be- 
haviour of the inter-electron terms in the Hamiltonian. 
Necessarily, such an approach is limited by the accuracy 
of the expansions employed. In addition, this increases 
the complexity of the computational problem. Heyl & 
Hernquist [35 in 1998 described both an analytical as 
well as a numerical approximation for evaluating the ef- 
fective inter-electronic potentials that was not only in- 
tuitive but also computationally less expensive. Their 
numerical technique extended the idea of adopting a ba- 
sis of functions for the directions both transverse and 
parallel to the field. They constructed the wave function 
in the axial direction with the help of harmonic oscilla- 
tor Hermite polynomials. This method was seen to yield 
accurate results for high magnetic fields j3z > 10 3 for 
hydrogen and helium. They employed a method of cal- 
culating the total energy of the system given an assumed 
set of wave functions and then proceeded to minimise this 
energy by varying the free parameters of the wave func- 
tions. This was seen to yield accurate results consistent 
with other work [BJ. One of the advantages over other 
methods [HI HO H3 El El (Ml 131] was observed to be 
the significantly lesser amount of computation involved. 
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More recently, Mori & Hailey [35] and Mori fe Ho [37] 
adopted a perturbative approach to treat the exchange 
terms and higher Landau states with success for finding 
upper bounds for the energies for the ground and first 
few excited states of helium and other mid-Z atoms in 
high magnetic field strengths. 

In the literature there are numerous studies outlining 
accurate estimates of the energies of different states of the 
hydrogen atom, however the number of investigations of 
the helium atom and helium-like species in the strong or 
intermediate magnetic field regime is rather small. The 
estimates of the energy levels of these species are only 
moderately accurate, and the computational expense is 
rather high. However, for most observable neutron stars 
and many white dwarf stars the magnetic field strengths 
lie in the intermediate field regime [BJ. In order to fa- 
cilitate a proper understanding of the spectra of neutron 
stars and white dwarf stars one must necessarily have 
more stringent bounds on the energy levels of atoms in 
the atmospheres of these compact objects in the interme- 
diate regime of magnetic field strengths. This is the aim 
of the current work. The work described herein extends 
previous work [351 and presents a numerical treatment 
of hydrogen and helium atoms in magnetic fields, yield- 
ing accurate results for the eigenvalues and eigenvectors 
of the first few low-lying states over a wide range of field 
strengths in the intermediate field regime. The calculated 
energy eigenvalues are seen to be improvements upon 
previous estimates. The procedures described herein do 
not make any assumptions of basis of functions and nei- 
ther are they restricted to the adiabatic approximation. 
The direct and exchange interactions of the electrons are 
computed using a novel method and do not rely upon 
any approximations. The method is readily extendable 
to many-electron systems and arbitrary field strengths. 
The overall method is also computationally straightfor- 
ward to implement. 
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It can be seen immediately that such a choice precludes 
the use of a basis of functions for describing the behaviour 
of the electron both parallel and perpendicular to the 
magnetic field. Thus, the time-independent Schrodinger 
equation in units of Bohr radii can be written 
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where m is the azimuthal quantum number. In defining 
Eq. (3) it has been assumed that the electron spin is 
anti-aligned with the magnetic field. Additionally, the 
energy parameter e, is defined as 
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Eryo is the Rydberg energy. The quantity a fts 1/137 is 
the fine structure constant. The parameter f3 is defined 
in the usual way as 
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where Bo is the critical field strength at which point the 
transition to the intense magnetic field regime occurs [6 . 
At this field strength the Larmor and the Bohr radii are 
equal in magnitude. 



II. THE HYDROGEN ATOM 

We shall begin with the Hamiltonian for the hydrogen 
atom in a magnetic field. The Hamiltonian for the single 
electron of the hydrogen atom in cylindrical coordinates 
is given by 
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where m e is the mass of the electron and the B the mag- 
netic field strength; the vector B is oriented along the 
positive z— axis. The remaining symbols have their usual 
meanings. It is of course implicitly assumed herein that 
the nucleus is infinitely massive. Let us assume a certain 
form for the wave function of the single electron 
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In general, the Larmor radius is related to the Bohr 
Radius via the equality an = as/3 -1 / 2 . Where an and 
as are the Larmor and Bohr radii respectively. Thus, 
beyond a value of (3 ~ 1 the transition to the intense 
magnetic field regime occurs and the interaction of the 
electron with the nucleus becomes progressively less dom- 
inant as (3 increases. The Larmor radius is given by 
an = \J eh/2m e while the Bohr radius is given by the 
expression ag — h/am e c. Based upon the above defini- 
tion of f3 in Eq. (5) it is convenient to classify the field 
strength [38] as low (/3 < 1CP 3 ), intermediate, also called 
strong (1CP 3 < /3 < 1), and intense or high (1 < (3 < oo). 

It can be seen that Eq. (3) is a linear second order par- 
tial differential equation, and in the present study it was 
solved numerically on a computer using finite-element 
techniques. For details on the numerical treatment of 
Eq. (3) see §4. 
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III. THE HELIUM ATOM 

For calculating the atomic structure of two-electron 
systems we adopt an iterative method, the so called 
self- consistent field method [27 which essentially solves 
the Hartree-Fock equations for the electrons. A short 
derivation of the key equations are given below assum- 
ing a single-configuration form for the atom as described 
above. 

It is to be noted that no restrictions are imposed upon 
the atomic electrons such as the commonly employed adi- 
abatic approximation [5J. 



Let us assume that the wave function of a given con- 
figuration of electrons is given by 



* = Ajy (V>i,V>2,V>3,— ,iPn-i,i/>n) (10) 

where, An is the anti-symmetrization operator. Thus, it 
can be seen that a single slater determinant is assumed to 
represent the atomic configuration of all electrons. The 
single particle wave functions are assumed to be of the 
same form as assumed for the case of the hydrogen atom 
in Eq. (2); explicitly, 



A. Derivation of the generalized Hartree-Fock 
equations in partial differential form 

Let us begin with the Hamiltonian of an ./V-electron 
atom split into one- and two-body terms 
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The first part of the Hamiltonian consisting of one- 
body interactions is given by the standard prescription 
(using polar cylindrical coordinates) 
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where i — 1, 2, N 



And the two-body term in the Hamiltonian is simply 
the Coulomb interaction between the i th and j th elec- 
trons. 
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*i = il>i(pi,Zi)e im + , Xi{'i) (11) 
where i labels the electrons for each of the N electrons. 
The single particle wave functions ipi(Pi, z i) are taken to 
be real functions. 

We shall look at only the fully spin-polarised states 
with each electron's spin being anti-aligncd with the mag- 
netic field. This was considered to be a reasonable as- 
sumption as the magnetic fields in white dwarf stars and 
neutron stars are high enough for such states to be the 
most tightly bound states [B]. 

Writing the generalized Hartree-Fock equations for de- 
termining the single particle wave functions tpi we have 



h{ri)tpi {r i )+^2[{i>j{r j )\w(r i ,r j )\<tp j (r j ))i) i (r i ) 
-{4>j(rj)\w {n ,rj)\ V>j (r,- ) } ip 3 (n)] 

= E^n) (12) 

where i = 1,2, 3, N 

Substituting the ansatz given in Eq. (11), the assumed 
individual electron wave functions, into Eq. (12) we ob- 
tain after rearranging some terms 
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where i, j — 1, 2, 3, N 



It is to be noted that the contribution due to electron 
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spin has been averaged out a priori. We have chosen 
to work in units of Bohr radii along with the definitions 
given below. Additionally, hereafter the exponential fac- 
tors with imcj) are to be interpreted with the appropriate 
sign depending upon whether they are written in the bra 
or in the corresponding ket; a minus sign for the former 
and a plus sign for the latter. This interpretation should 
presently be applied to Eq. (13). 

The Bohr radius for an atom of nuclear charge Z is 
given by a^/Z, with as as defined earlier. The magnetic 
field strength parameter is defined as flz = B/Bz = 
(3/Z 2 . The reference magnetic field strength for nuclear 
charge Z, at which the transition to the intense mag- 
netic field regime occurs is given by Bz = Z 2 Bq with 
B as defined in Eq. (6). Finally, the energy parame- 
ter is defined as £j = Ei/(Z 2 E OQ ), with Eoo as defined 
in Eq. (4). The above written Eq. (13) represents the 
N-coupled Hartree-Fock equations in partial differential 
form for an iV-clcctron system with nuclear charge Z. 
The system of equations is solved iterative ly; see §4 for 
numerical details. 

Of key concern in the computation of the eigenvalues 
and eigenvectors is the determination of the direct and 
exchange interactions between the electrons. In the cur- 
rent study these have been dealt with in a manner rather 
different from earlier treatises. In short, these contribu- 
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Which then immediately yields 

V.H>d = -4^(p^)| 2 (16) 

The RHS of Eq. (16) is the square of the j th electron's 
wave function evaluated using the coordinates of the i th 
electron. Noting that Eq. (16) is simply the Laplace 
equation, it is observed that it is numerically tractable 
and solved using the appropriate boundary conditions 
to yield the potential $£> which is due to the direct in- 
teraction between the i th and j th electrons; see §4 for 
details on the numerical methods employed. It is to be 
noted however, that in contrast to previous work, the 
problem in Eq. (16) is somewhat simpler despite hav- 
ing to solve a partial differential equation as one does 
not have to find approximate expressions for the mixing 
terms arising from the interaction between different elec- 
tronic states and the Coulomb potential. The reader is 
referred to [6j [29] and references therein for the different 



tions are essentially extra potentials that add to the exist- 
ing single particle Hamiltonian, except they are coupled 
through the wave functions weighting them in a given 
equation in Eq. (13). The procedure employed for eval- 
uating these potentials is described below. 



B. The direct interaction 

The treatment of the direct and exchange potentials 
considered here and in the subsequent section is based 
upon rigorously determining these potentials by solving 
their corresponding elliptic partial differential equations. 
This is carried out in this study by employing the square 
of the individual electrons' momentum operators — ifiVj. 
The details of method developed in this study are given 
below. 

Let us first examine the integral representing the direct 
interaction between the electrons. 
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Let us act on both sides of Eq. (14) with the operator 
p 2 = — /i 2 Vf , we then obtain 



-4n6 3 (r i -f' j )\i }j (p j ,z j )e im ^) (15) 
I 

approximation methods employed in obtaining estimates 
for Eq. (14). 

We now turn our attention to the other two-body term 
in Eq. (13); the exchange interaction. 

C. The exchange interaction 

We shall follow the same methodology as in our treat- 
ment of the direct interaction term. Let us re-write the 
term in Eq. (13) that relates to the exchange interaction 
between the i th and j th electrons 



<S> E = (^( PjiZj ) e ^\—--\Mp J ,z J )e^) (17) 

I M T j I 

Again, as in our previous treatment, let us act on both 
sides of Eq. (17) with the operator Vf , this time drop- 
ping the redundant factor of — h 2 , to obtain 
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Upon carrying out the integral we get 
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It is to be noted that by the definition in Eq. (11) ip* — ip for the spatial part of the individual electron wave 
functions. At this stage, let us make the ansatz that 



^E = a E (p i ,z i )e i ^- m ^ i 
Let us act on both sides of Eq. (20) with the Laplacian operator, Vf , to obtain 
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It is then a straight-forward matter upon comparing Eq. (19) with Eq. (21) to immediately see that 
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The elliptical partial differential equation, Eq. (22) is solved numerically and we thus obtain an estimate for the 
function a E (pi,Zi) for each of the pair-wise interactions among the TV electrons. Knowing a E (pt,Zi) we can then 
obtain according to Eq. (20). Once <&e and <!>£> have been obtained we can substitute them into Eq. (13) for the 
potentials due to the direct and exchange interactions respectively to obtain 
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Eq. (24) is the final form for the generalized Hartree- 
Fock Equations and for two-electron systems we have two 
equations, however for arbitrary nuclear charge Z. The 
energy of the Hartree-Fock state is given by Eq. (25). 



This completes our derivation of the Hartree-Fock 
equations for atoms in magnetic fields of arbitrary 
strength. The following section delineates the numeri- 
cal procedures employed in the calculation of the energy 
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eigenvalues and eigenfunctions. Thereafter, results are 
presented and a discussion follows. 



IV. NUMERICAL DETAILS 

The eigenvalue problem for the hydrogen atom in Eq. 
(3) is solved by discretising the equation and solving the 
resultant algebraic eigenvalue problem. The discretisa- 
tion is done using the finite-element method (FEM) [39] . 
The generalized eigenvalue problem is then solved using 
a sparse matrix generalized eigensystem solver [40 . This 
method was found to yield accurate results for the energy 
eigenvalues for the first few cigcnstatcs with different az- 
imuthal quantum number, m;. Runs were carried out for 
different values of the magnetic field strength parameter 
(3 in the range 10 -2 < j3 < 10. For testing additional con- 
vergence for every run we employed six different levels of 
mesh refinement, ranging from coarse to sufficiently fine 
mesh. The fine mesh calculations took up to two days of 
computing time on AMD Opteron® 844 1.8 GHz pro- 
cessors. 

For solutions to the helium atom in strong magnetic 
fields an atomic structure software was developed as a 
part of this study for the purpose of calculating the en- 
ergies of different states of multi-electron atoms. The 
program takes as its input, the number of electrons in 
the atom n e , the nuclear charge Z, and the magnetic 
field strength parameter (3 and then proceeds to compute 
systematically the eigenvalues and eigenfunctions of the 
coupled system of equations in Eq (24) according to the 
iterative procedures described below in brief. 

The Eqs. (16), (22) and (24) are solved in a three 
step process. First an initial estimate is obtained for 
the eigenvectors by solving Eq. (24) without the contri- 
butions due the interaction between the electrons. The 
wave functions are then used as starting points for further 
improvements. The second step involves obtaining esti- 
mates for the potentials due to the direct and exchange 
interactions amongst the electrons, vis-a-vis, the ellipti- 
cal partial differential Eqs. (16) and (22) are solved using 
the estimates for the wave functions obtained in the pre- 
vious step. These estimates are then used to solve for 
better estimates of the eigenfunctions along with the rel- 
evant eigenvalues in Eq. (24). The last two steps are iter- 
ated in the order described above to obtain progressively 
better estimates for the eigenvalues and eigenvectors with 
each iteration. It was observed during our runs that fast 
convergence was achieved; within the first few iterations. 
A convergence criterion was employed wherein the differ- 
ence between the solutions for two consecutive iterations 
was tested. If the difference was smaller than the desired 
tolerance for convergence then it was considered to be a 
solution of the iterative scheme. Typically, a tolerance 
on the order of 10 -6 was employed. Thereafter, the to- 
tal energy of the Hartree-Fock state under consideration 
is reported according to Eq. (25). Again as in the case 
of the hydrogen atom runs were carried out for different 



values of (3 with five different levels of mesh refinement 
in each case for additional convergence testing. 

Once again, the generalized system of eigenvalue equa- 
tions were solved using appropriate FEM discretisation 
yielding accurate results for the eigenvalues and eigen- 
vectors for the first few low-lying states of helium. For 
details on the finite-clement method the reader is referred 
to Ref. [3S]. It was found that domain compactification 
was computationally expensive in terms of memory and 
computation time for obtaining sufficiently accurate re- 
sults. On the other hand, limiting the domain size in the 
two orthogonal directions (parallel and perpendicular to 
the magnetic field) to several Bohr radii (~ 20) was ob- 
served to give accurate results for the range of magnetic 
field strengths considered in this study. Concordantly the 
computational expense was many times less in compar- 
ison to the former case. Numerical errors arising from 
truncating the domain were not significant {i.e. they 
were smaller than the discretization errors). 

V. RESULTS & DISCUSSION 

The results from the calculations carried out for hydro- 
gen and helium atoms in strong magnetic fields are pre- 
sented below. Only the results pertaining to the ground 
and first two excited states of these atoms in strong 
magnetic fields are given. It is to be mentioned that 
the present study focuses on calculations for magnetic 
field strength within the so called strong or intermediate 
regime (10 -2 < /? < 10). 

A. The Hydrogen Atom 

For magnetic fields in the range 10~ 2 < (3 < 10 the 
time-independent Schrodinger equation given in Eq.(3) 
was solved for different values of the magnetic field 
strength parameter. However only the results obtained 
for the most tightly bound states with azimuthal quan- 
tum number, m = 0, —1 and —2 are presented here. The 
binding energies are reported in Rydberg units. (The 
states with m > become unbound with increasing mag- 
netic field strength [6].) 

The variation in the binding energy for the state cor- 
responding tom = 0,7r = +lis shown in Figure [T] The 
quantity ir indicates parity with respect to the z— axis. 
The data points arc eigenvalues obtained from the numer- 
ical solution of Eq. (3). The energy eigenvalues are re- 
ported as the values corresponding to infinitely fine mesh 
sizes, or in other words when the average area of the fi- 
nite elements approaches zero. A short discussion of this 
estimation procedure is given later. As can be seen in the 
figure, the electron becomes more and more bound as the 
magnetic field strength increases. The line through the 
data represents a fit to the data. A rational function 
was used to model the data in this regime using a ro- 
bust Levenberg-Marquardt method jH]. The values for 
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TABLE I: Table listing the coefficients of the different rational 
functions employed for fitting the three states of hydrogen 
discussed. The maximum fractional error of the eigenvalue 
relative to the fit from f3 = 10~ 4 to f3 = 10 3 follows the the 
list of coefficients. 

State Coefficients State Coefficients 



a = 3416640.9 





a 


3.0849526 




ai = 


12197155 




a = 


6.6571666 




a 2 = 


4879198.1 


ls 


a 2 = 


—0.14904847 


ls 


a-i = 


425959.45 


a.3 = 


0.43176735 


0,4 — 


119957.96 




bo = 


3.0792681 
9 x 10~ 3 




bo = 
bi = 


3418004 
5216732.1 
4 x 10" 4 




do = 


0.58201924 




Oo = 


0.038324719 




ai = 


13.769289 




Oi = 


2.7226023 




a 2 = 


17.780336 




a 2 = 


6.9791245 


2p_i 


a.3 = 


-0.52325377 


3d-2 


a3 


-0.55680026 


0,4 = 


0.74135964 


a ; = 


0.4177192 




bo = 


2.3266637 




bo = 


0.34444067 




bi = 


18.424245 
3 x IO -3 




bi = 


6.7184551 
1 x 10 -2 



the coefficients of the interpolating function are given in 
Table i 

These rational functions accurately model the data in 
a range of magnetic fields from (3 = 10 -4 to (3 = 10 3 
(twice the range indicated in Figure [lj and could poten- 
tially be used directly in atmosphere models for neutron 
stars and white dwarf stars. Having such accurate ana- 
lytical forms for the energies of atoms in strong magnetic 
fields obviates the need for performing laborious calcula- 
tions for the energy estimates making atmosphere models 
computationally less intensive. In addition, it becomes 
possible to analyze observed spectra of neutron and white 
dwarf stars with relative ease at arbitrary field strengths 
within the intermediate field regime. The fits to the data 
are rational functions obtained as described above. The 
functions have the form 



f(x) = l \ (26) 

where x = ln(l + /3). The coefficients and the maximal 
fitting errors over the entire range (3 = 1CP 4 to (3 = 10 3 
are given in Table [i] The data in the range [3 — 10 2 to 
(3 = 10 are results from calculations of the present study 
while, to construct the fits for larger and smaller values of 
(3 than the indicated range we used the results of Rcf . [6 1 . 

As can be clearly seen in the figure the binding en- 
ergy of the state of hydrogen with quantum numbers 
m = 0, 7r = +1 increases monotonically with increasing 
magnetic field strength. The dashed line in the graph 
represents a first order perturbation theory calculation. 
The purpose being to illustrate the fact that perturba- 
tion theory breaks down with increasing the magnetic 
field strength. It is evident upon inspection that the 
breakdown of perturbation theory occurs within the so 
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FIG. 1: Figure showing the variation in the binding energy 
of the ground state of hydrogen m = 0, n = +1 with the 
magnetic field strength parameter in the range 10" 2 < (3 < 
10. The energy is reported in units of Rydberg energy. The 
data points are results of the numerical calculation while the 
line represents a fit to the data. The dashed line represents 
the results from first order perturbation theory. 



called intermediate field regime. Thus arises the need for 
accurate data for the structure of atoms in this regime of 
magnetic field strength 10 -2 < f3 < 10. 

Figure [2] shows variation in the binding energy the 
three most tightly bound states of the hydrogen atom 
with magnetic field strength, [3. Explicitly, in order of de- 
creasing binding energy, the three states m — 0; 7r = +1, 
m = — l;7r = +1 and m = — 2;7r = +1, arc plotted in 
Figure [2] Again the lines through the data points rep- 
resent fits to the data. Every data point in these figures 
was obtained as an estimate corresponding to limit of the 
finite element size going to zero; see later. It can be seen 
in Figure [2]that the binding energy of the different states 
increases dramatically after (3 ~ 1. Thereafter, the bind- 
ing energy increases at an increasing rate with increasing 
magnetic field strength. 

It is to to be noted that, in all cases, the rational func- 
tions were so chosen as to reflect the fact that eventually 
at large values of the magnetic field strength parameter 
(3 the energy would be proportional to In 2 j3 [8] . For the 
sake of comparison with previously determined accurate 
results, the fit to the calculated data for the ground state 
of hydrogen corressponding to m = 0, 7t = +1 is plotted 
against the calculated data of the researchers in Ref [6] in 
Figure [3] The line through their data is a spline fit. As 
can be seen in the figure, the two lines are nearly com- 
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FIG. 2: Figure showing the variation in the binding ener- 
gies of three lowest m- states of hydrogen corresponding to 
m = 0, — 1, —2, with n — +1 with the magnetic field strength 
parameter in the range 10 -2 < f3 < 10. The energy is re- 
ported in units of Rydberg energy. The data points are re- 
sults of the numerical calcuation while the lines represent fits 
to the data. 



pletely coincidental. The differences were seen to be in 
the fifth or sixth decimal place between the results of the 
present study and those of Ruder et al [5J; or in other 
words the results were considered to be of comparable 
accuracy. For details of the energy eigenvalues obtained 
in this study the reader is referred to Table |III| in the 
appendix to this paper. 

Figure [4] shows how the eigenvalues obtained from the 
solution of Eq. (3) depend upon the mesh size. As was 
mentioned earlier, the calculations were carried out on 
a finite domain of several Bohr radii in each of the two 
directions both parallel and perpendicular to the field. 
While keeping the domain fixed, the number of finite 
elements constituting the mesh was varied. Runs were 
performed on each mesh size for every value of magnetic 
field strength considered in the study for each of the three 
states of hydrogen, m — 0, —1 and —2. The eigenvalues 
thus obtained were plotted against the average area per 
finite element in the mesh corresponding to each value 
of the magnetic field strength parameter (3. Extrapola- 
tion of the data to the limit of zero mesh size in each 
case yielded the values that would correspond to an in- 
finitely fine mesh. These values were reported as the cal- 
culated data points in the preceding figures. The extrap- 
olation was carried out by employing rational functions; 
the reader is referred to Ref. [H] for details regarding the 



FIG. 3: Figure showing the variation in the binding energies 
of the state of hydrogen corresponding to m = with -n = 
+ 1 with the magnetic field strength parameter in the range 
10~ 2 < j3 < 10. Results from both the present study as well as 
those from Ref [6] are plotted herein. The energy is reported 
in units of Rydberg energy. The line through the data from 
Ref. [6] is a spline fit. With regard the current study, the data 
points are results of the numerical calcuation while the line 
represents a fit to the data. 



method employed. The average accuracy of the estimate 
of the asymptotic value was determined to be on the or- 
der of 3 x 1CP 6 Rydbergs. Such estimation errors were 
too small to be shown on the plots. The line through the 
data points is merely a guide to the eye. 

The wave functions for the most tightly bound state 
of hydrogen to = 0, 7r = +1 are plotted in Figure [5] 
for four different values of the magnetic field strength 
(3. The plots represent a slice through the three dimen- 
sional atom in the upper half of the r — z— plane. The 
complete representation in three dimensions is the figure 
of revolution about the z— axis and simultaneously re- 
flected about the xy— plane. The length units are Bohr 
radii of the hydrogen atom and the a:— axis represents 
the direction perpendicular to the magnetic field direc- 
tion while the y— axis represents the direction parallel to 
the magnetic field, i.e. the z— direction in three dimen- 
sional cylindrical polar coordinates. The purpose here is 
to illustrate the dramatic change that occurs near /3 rs 1 
when the electron becomes tightly bound and the bind- 
ing energy increases dramatically with increasing [3. It is 
immediately evident upon inspection that the spherical 
symmetry of the atom is clearly broken as we approach 
higher and higher magnetic field strengths. 



9 



"l l [ l [ [ l l | 




I I 1 1 I 1 I I I 1 I 

io- 3 io- 2 io- 1 10° 

Area Per Finite-Element [(Bohr Radii) 2 ] 



FIG. 4: Figure showing a plot of the eigenvalues calculated 
as a function of the area per finite element in the mesh for 
a calculation for the ground state of hydrogen at /3=10 -2 . 
With increasingly finer meshes for a sufficient domain size the 
eigenvalue obtained approaches an asymptote corresponding 
to the limit of infinitely fine mesh. The parameter describing 
the mesh size is the average area per finite element in the 
mesh 



B. The Helium Atom 

For magnetic fields in the range 10" 2 < 0z < 10 we 
solved the set of coupled Hartree-Fock equations given in 
Eq.(24) for different values of the magnetic field strength 
parameter. We shall however be showing the results ob- 
tained for the states corressponding to M = 0, —1 and —2 
corresponding to the three most tightly bound states of 
helium. The two electrons are taken to posess azimuthal 
quantum numbers given by m\ and m 2 respectively with 
M = m\ + mi. We consider here in the study only the 
fully spin-polarised states, thus each of the electrons has 
its spin anti-aligned with the magnetic field as assumed 
in the derivation of Eq. (24). These states have been seen 
to be the most tightly bound states of helium in strong 
magnetic fields [B]. It is to be noted that the binding 
energies are reported in units of Rydberg energy in the 
Coulomb potential of nuclear charge Ze (alternatively 
called atomic units or simply a.u.) 

At the very outset, as a preliminary test of the Hartree- 
Fock atomic structure software developed as part of this 
study, we carried out calculations for the binding energy 
of the 1 So singlet state of neutral helium without any 
magnetic field (B = case). However, in this case the 
wave function of the configuration of electrons is com- 
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FIG. 5: Figure showing a plot of the change in the real part of 
the wave functions, i.e. ip(p, z), of the ground state of the hy- 
drogen atom with varying magnetic field strength. The wave 
functions show the ground state of hydrogen with quantum 
numbers m = 0, s z = —1/2 and n = +1 in magnetic fields 
of strength corresponding to (a) f3 = 10 -2 , (b) (3 = 10 _1 , (c) 
f3 = 1 and (d) [3 = 10. It can be clearly seen that as the 
magnetic field strength increases so does the binding energy 
as indicated by the clear shrinking of the wave function with 
increasing (3. These plots are not to be confused with electron 
density. 



pletely symmetrized with respect to the spatial part of 
the total wave function while the spins of the two elec- 
trons are anti-parallel to each other. With these changes 
Eq. (24) was solved using the numerical procedures out- 
lined in §4 and the energy reported according to Eq. (25) . 
The runs were carried out for different mesh sizes over the 
domain and the final eigenvalue is reported after extrap- 
olating to the limit of infinitely fine mesh as described 
above. The eigenvalue obtained using this procedure was 
1.4499 au., while the most accurately determined result 
via numerical techniques thus far is 1.4519 au The 
difference is about 2 x 10~ 3 au. This was considered 
to be sufficiently accurate given the fact the calculations 
carried out as part of this study were single-configuation 
calculations while the result from Rcf. [42 is essentially a 
multi-configuration calculation which is a computation- 
ally more intensive method. Thus, the atomic structure 
software developed for this study was considered to be 
sufficiently accurate for the purposes of the current study. 

The variations in the binding energies for the most 
tightly bound states corresponding to M = —1,-2 and 
with S z — —1 are shown in Figure [6] The data points 
are eigenvalues obtained from the numerical solution of 
Eq. (24) according to the numerical procedures described 
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in §4. The energy eigenvalues are once again reported as 
the values corresponding to infinitely fine mesh sizes, es- 
timated according to the discussion above; see Figure [4] 
and discussion thereof. As in the case of hydrogen, runs 
were performed on each mesh size for every value of mag- 
netic field strength considered in the study for each of the 
three states of helium, M = 0,-1 and —2. These eigen- 
values were extrapolated using rational functions [41 to 
the limit of zero mesh size. Eigenvalues obtained by this 
extrapolation procedure are reported as the calculated 
data points in Figure [6] The average accuracy of the es- 
timate of the asymptotic value was determined to be on 
the order of 2 x 10~ 5 Rydbergs. Again, for details re- 
garding the extrapolation method, the reader is referred 
to Ref. [H]. As can be seen in Figure [6j the binding 
energy of the state increases with increasing magnetic 
field strength (3z- The line through the data represents 
a fit to the data. In addition, a good measure of the 
discretization error was considered to be the difference 
between the computed eigenvalues for the most finely re- 
fined mesh size employed and the extrapolated result for 
the mesh size tending to zero. This error associated with 
the discretization of the problem is reported in the num- 
ber appearing in the parentheses in Tables IV V] and VI 



in the appendix to this paper; the number therein cor- 
responds to the absolute error in the fifth decimal place. 
Concordantly, the error bars were too small to be shown 
on the plots. 

Again a rational function was used to model the data in 
this regime using a robust Levenberg-Marquardt method 
|41j . The values of the coefficients of the interpolating 
functions are given in Table |TTJ It is to be noted at this 
juncture that the rational functions were so chosen as 
to reflect the fact that eventually, for large values of the 
magnetic field strength parameter (3z, the binding ener- 
gies are proportional to In 2 fiz- With this in mind, the 
aim was to find an accurate analytical function that de- 
scribed the variation in binding energy over a wide range 
of (3z such that it employed the least number of parame- 
ters. Inspection of Table [TT| reveals that this was accom- 
plished with five free parameters for modeling the data in 
the range of magnetic field strength 10 -2 < (3z < 10. Po- 
tentially, as was noted earlier, these analytical functions 
can be employed with relative ease in atmosphere models 
directly, thus circumventing the need for both laborious 
calculations of the energies via HF methods and simul- 
taneously avoiding the need for spline interpolations of 
tabulated data of binding energies. The maximum errors 
in the fits are also provided in Table [TTJ 

As can be clearly seen Figure [6] the binding energies 
increase monotonically with increasing magnetic field 
strength. Beyond (3z ~ 2.5 x 10 _1 which corresponds 
to (3 w 1 the departure from spherical symmetry of the 
atom becomes more pronounced as the magnetic field 
strength increases; see Figure 10 This is indicated by 




the steep rise in the curves signifying the increasing rate 
at which the binding energies increase with increasing 
magnetic field strength. 



FIG. 6: Figure showing the variation in the binding energy of 
the three most tightly bound states of helium with the mag- 
netic field strength parameter in the range 10" 2 < (i z < 10. 
The energy is reported in units of Rydberg energy for nu- 
clear charge Ze. The data points are results of the numerical 
calculation while the lines represent fits to the data. 



For the sake of comparison with previously determined 
results, the fits to the calculated data for the three 
most tightly bound states of helium corresponding to 
M = —1,-2 and with S z = — 1 are plotted against 
the calculated data of the researchers in Ref. [6] and 
Refs. [321 [34] in Figures [7] [8] and [9j The lines through 
the data of the other researchers are spline fits. It can be 
seen that the rational functions used to fit the calculated 
eigenvalues for the binding energies of the three states of 
the helium atom are consistent with previous work. In 
particular they are seen here to be improvements upon 
the estimates of Ruder et al [5] over the entire range of j3z 
considered in this study. Jones et al on the other hand 
had employed a Monte-Carlo approach [55] for solving 
the HF equations. They assumed a large number of ba- 
sis functions with variable parameters that could be fine 
tuned within the framework of Monte-Carlo simulations 
to arrive at upper bounds for the energies. Though this 
method is effective, as can be seen from their estimates 
shown in Figures [7] [8] and [9] it is computationally de- 
manding and it restricts the wave functions of the elec- 
trons to be expressed using a finite number of basis func- 
tions. The method described in the current study does 
not impose such a condition and thus the wave functions 
that are determined are in effect superpositions of a very 
large number of such basis functions and arise naturally 
from the solution itself. However, it is to be mentioned 
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that for high magnetic field strength parameter [3z ~ 10 
for the state M = -1 and (3 Z ~ 20 for the states M = -2 
and 0, the eigenvalues were seen to drop below those of 
Jones et al's estimates; see Figures [7] [8] and [9] It is 
to be noted that this effect was due to insufficient com- 
puter memory to carry out calculations with the desired 
mesh refinement at high magnetic field strengths where 
the electrons are tightly bound to the nucleus and their 



wave functions shrink closer to the nucleus; see Figure 10 



With sufficient computer memory it is expected that the 
results presented here can be extended to higher mag- 
netic field strengths. 

In Figure [7] it can be seen that the estimates of the 
present study are improvements upon the estimates of 
Ruder et al |§J. The long-dashed line representing the 
data of the researchers in Ref. [6] has a kink at about 
fiz = 1, this is because they used cylindrical rather than 
spherical coordinates at this stage for solving the eigen- 
value problem. The eigenvalues obtained are seen to be 
better estimates than those of Jones et al [351 [M] in the 
range 10~ 2 < /3z < 2 x 10 _1 by about 7% maximum 
and 4% minimum. However above that field strength, 
the estimates of the current study fall below those of 
Refs. [521 154] by about 2% maximum. Since the present 
study is a numerically less complicated technique, thus it 
was expected that the Monte-Carlo simulations of Jones 
et al [32 , 34 would arrive at a result slightly better than 
those of the current study given the limitations of com- 
puting memory at higher magnetic field strengths. How- 
ever, one of the aims of the current work to was facilitate 
a numerically straightforward and computationally less 
intensive method for computing the energies of atoms in 
strong magnetic fields and in that regard the results ob- 
tained herein were considered to be good. 

In Figures [8] and [9j it can be seen that estimates of the 
binding energies of the first and second excited states of 
the helium atom in strong magnetic fields are consistent 
with the findings of other researchers [5J [35J 151] . The 
range of improvements were between 11% maximum to 
about 1% minimum over the entire range 10 -2 < (3z < 10 
relative to the eigenvalues obtained by the the researchers 
in Refs. [32, 34 . Towards higher magnetic field strengths 
it was observed that the improvements in the estimates 
tended to drop; this is due to the fact that at such high 
magnetic field strengths for obtaining more accurate es- 
timates of the energies, greater computer memory was 
required to accommodate for finer meshes and this was 
not possible within the framework of the current project. 
On the other hand, it can be seen that the results of the 
current work presented in Figures [7] [8] and [9] are seen 
to be on average 8% improvements relative to the re- 
sults obtained by Ruder et al [6] ; the reader is referred to 
Tables [TV] |V| and | VI[ provided in the appendix to this pa- 
per, for exact eigenvalues of the three most tightly bound 
states of helium calculated in this study. A comparison 
with the results of other researchers can also be found 
therein. Once again the kinks in the curves in Figures [8] 
and [9] of the results of Ruder et al [H] at around [3z ~ 1 
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FIG. 7: Figure showing the variation in the binding energies 
of the state of helium corresponding to M = — 1 with 7r z = 
+ 1 with the magnetic field strength parameter in the range 
10~ 2 < Pz < 10. Results from both the present study as 
well as those from Ref [6] and Jones et al [321 134] are plotted 
herein. The energy is reported in units of Rydberg energy 
for nuclear charge Ze. The data points are results of the 
numerical calcuation while the lines represent fits to the data. 



are due to the fact that they changed from a spherical to 
a cylindrical representation of the wave functions. 

Finally, the wave functions for one of the electrons in 
the most tightly bound state determined from the calcu- 
lations are plotted in Figure [10] for four different values 
of the magnetic field strength f3z- For the purpose of 
illustrating the dramatic change in the structure of the 
atom with increasing magnetic field strength we chose 
to show the electron in the state with quantum numbers 
m = — l,s z = — 1/2, 7r z = +1, which is equivalent to 
2p_i orbital in the low field limit. The length units are 
Bohr radii of the helium atom, i.e. clq/Z and the x— axis 
represents the direction perpendicular to the magnetic 
field direction while the y— axis represents the direction 
parallel to the magnetic field, i.e. the z— direction in 
three dimensional cylindrical coordinates. It can be seen 
that the electron becomes more and more tightly bound 
and the binding energy consequently increases dramati- 
cally with increasing f3 z . It can be seen that comparing 
Figure 10 'd) with Figures 



10'a-c) that the electron wave 



function shrinks considerably over the range of magnetic 
field strengths considered in the current study. 
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FIG. 8: Figure showing the variation in the binding energies 
of the state of helium corresponding to M = —2 with 7r z = 
+1 with the magnetic field strength parameter in the range 
1CP 2 < (iz < 10. Results from both the present study as 
well as those from Ref [6] and Jones et al [321 134] are plotted 
herein. The energy is reported in units of Rydberg energy 
for nuclear charge Ze. The data points are results of the 
numerical calcuation while the lines represent fits to the data. 



FIG. 9: Figure showing the variation in the binding energies 
of the state of helium corresponding to M = with n z = 
— 1 with the magnetic field strength parameter in the range 
10~ 2 < Pz < 10. Results from both the present study as 
well as those from Ref [6] and Jones et al [321 134] are plotted 
herein. The energy is reported in units of Rydberg energy 
for nuclear charge Ze. The data points are results of the 
numerical calcuation while the lines represent fits to the data. 



TABLE II: Table listing the coefficients of the different ratio- 
nal functions for fitting the three states of helium discussed. 
The maximum fractional error of the eigenvalue relative to 
the fit from (3 Z = 10~ 2 to (3 Z = 10 1 follows the the list of 
coefficients in each case. 



State 




Coefficients 


State 




Coefficients 




do 


= 0.19161951 




ao 


= 0.56302392 




(>i 


= 1.8523924 




en 


= 2.5417308 


M = -1 


a 2 


= 1.5767666 


M = -2 


0,2 


= 1.2729579 


lso; 2p_i 




= 0.10192465 


lso; 3d_2 


az 


= 0.1570202 




6 


= 0.1655611 




bo 


= 0.4756007 






1 x 10 -2 






6 x 10 -3 





ao = 


0.15433115 




ai = 


1.5812466 


M = 


a 2 = 


1.3219681 


lso;2p 


a-3 = 


0.034875245 




6o = 


0.13137227 






9 x 10~ 3 



VI. CONCLUSIONS 

The work described here was motivated by the need to 
have accurately determined values for the upper bounds 
for the energy levels of atoms in strong magnetic fields. 
As was discussed earlier, this need has arisen due to the 



presence of strong magnetic fields in neutron stars and 
white dwarf stars. The most commonly present atoms in 
the atmospheres of these compact objects, hydrogen and 
helium were studied here with the intention of obtaining 
accurate estimates of the energy levels of the first few 
low lying states in strong magnetic fields. We described 
a method adopting a physically motivated approach gov- 
erned by the inherent symmetries of the problem. We 
simultaneously circumvented the need for adopting a def- 
inite basis of functions to describe the wave functions of 
the electrons in either of the directions, parallel and per- 
pendicular to the magnetic field. The approach is un- 
restrictive with regard to the wave function; it has the 
distinct advantage over methods that require a basis of 
functions to describe the wave functions because in nu- 
merical solutions one can only have a finite number of 
such functions. The wave functions determined in the 
present study came about naturally from the symmetries 
of the problem and are thus in effect superpositions of a 
large number of basis functions. 

Such an approach however resulted in elliptical partial 
differential equations for the electrons that were subse- 
quently solved using finite clement techniques. It is to be 
noted that the novel method adopted for determining the 
direct and exchange interactions between the electrons in 
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1 (a) [3 Z = 10 



(b) f3 z = 10 _1 ,m 2 = -1 




10 10 

Distance across field [Bohr Radii] 



FIG. 10: Figure showing a plot of the real part of the wave 
functions, i.e. lp(p,z), of one of the electrons comprising the 
ground state of the helium atom with changing magnetic field 
strength. The wave functions show the outer electron in the 
most tightly bound state of helium with quantum numbers 
77i2 = — 1, s z = — 1/2, 7r = +1 in magnetic fields of strength 
corresponding to (a) f3z = 10" 2 , (b) f3z = 10" 1 , (c) fiz = 1 
and (d) f3z = 10. It can be clearly seen that as the magnetic 
field strength increases so does the binding energy as indicated 
by the clear shrinking of the wave function with increasing j3z- 
These plots are not to be confused with electron density. 



the helium atom is also exact in the sense that it does 
not rely upon any ab initio assumptions to approximate 



the integrals. These interaction potentials are solved in 
a natural manner by solving the elliptical partial differ- 
ential equations, Eqs. (16) and (22). The eigenvalues 
found in the range of the magnetic field strength param- 
eter 10 -2 < f3, f3z < 10 considered in this study were seen 
to be consistent with previous findings [6j (321 IS] • Ratio- 
nal functions were also used to find sufficiently accurate 
interpolating functions for the binding energies of vari- 
ous states of both the hydrogen and helium atoms in the 
range of magnetic fields considered herein. These were 
seen to be accuarate to (an average for all six fits) within 
0.8%. Potentially such interpolating functions could be 
used in atmosphere models of neutron stars and white 
dwarf stars thus obviating the need for involved and la- 
borious calculations of the same. 



Thus, the current work describes an unrestricted and 
computationally less intensive method for calculating the 
energy levels of atoms in strong magnetic fields. There 
are in essence three directions in which the current work 
could be extended. First, the current work can be readily 
extended to higher magnetic field strengths by changing 
the domain of calculations appropriately to incorporate 
the fact that the electrons become more and more bound. 
Simultaneously, the calculations and the software devel- 
oped as a part of this study are readily extendable to 
systems with more than two electrons, i.e., Li, Be, B, C, 
O etc. Finally, the procedures implemented herein can 
also be extended towards a multi-configuration frame- 
work [43 . In essence, the calculations employed herein 
are for a single configuration of the electrons. Thus hav- 
ing found the energies of the different configurations, it 
is possible to extend the theory to incorporate a multi- 
configuration approach which is likely to improve the re- 
sults already obtained here. 
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APPENDIX A: TABLES OF BINDING ENERGIES 



TABLE III: Table listing the binding energies of the three most tightly bound states of hydrogen. Energies are in units of 
Rydbergs. The states are labeled according to their counterparts in the low-field limit. Results from Ref. [6] are also provided 
for comparison. The number in the parentheses is the absolute error at the fifth decimal place; this is determined as the 
difference between the computed eigenvalues for the most finely refined mesh size employed and the extrapolated result for the 
mesh size tending to zero. 
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TABLE IV: Table listing the binding energies of the most tightly bound state of helium in moderate to large magnetic fields; 
M = —1,5*2 = — l,7r z = +1. Energies are in units of Rydberg energies in the Coulomb potential of nuclear charge Ze, where 
Z = 2 for helium. The results from the current work can be compared readily with previous work by the researchers in 
Refs. 6 31, 32, 34 . The number in the parentheses is the absolute error at the fifth decimal place; this is determined as the 
difference between the computed eigenvalues for the most finely refined mesh size employed and the extrapolated result for the 
mesh size tending to zero. 
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TABLE V: Table listing the binding energies of the low-lying state of helium; M = —2,S Z = — l,n z — +1. Energies are in 
units of Rydberg energies in the Coulomb potential of nuclear charge Ze, where Z = 2 for helium. The results from the current 
work can be compared readily with previous work by the researchers in Refs. 6, 32, 34 . The number in the parentheses is the 
absolute error at the fifth decimal place; this is determined as the difference between the computed eigenvalues for the most 
finely refined mesh size employed and the extrapolated result for the mesh size tending to zero. 
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TABLE VI: Table listing the binding energies of the low-lying state of helium; M — 0, S z — — 1, n z — — 1. Energies are in units 
of Rydberg energies in the Coulomb potential of nuclear charge Ze, where Z = 2 for helium. The results from the current work 
can be compared readily with previous work by the researchers in Refs. 6, 31, 32, 34 . The number in the parentheses is the 
absolute error at the fifth decimal place; this is determined as the difference between the computed eigenvalues for the most 
finely refined mesh size employed and the extrapolated result for the mesh size tending to zero. 
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